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ABSTRACT 


This report presents a comparison of four different estimation techniques 
applied to the problem of continuously estimating the parameters of a 
sinusoidal Global Positioning System (GPS) signal, observed in the presence of 
additive noise, and under extremely high-dynamic conditions* Frequency 
estimates are emphasized, although phase and/or frequency rate are also 
estimated by some of the algorithms* These parameters are related to the 
velocity, position, and acceleration of the maneuvering transmitter. Estimated 
performance at low carrier-to-noise ratios and high dynamics is investigated 
for the purpose of determining the useful operating range of an approximate 
maximum likelihood (ML) estimator, an extended Kalman filter (EKF), a 
cross-product automatic frequency control (CPAFC) loop, and a digital 
phase-locked loop (PLL). Numerical simulations are used to evaluate 
performance while tracking a common trajectory exhibiting high dynamics. 
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I. INTRODUCTION 


The problem of estimating the parameters of a single-frequency tone in 
the presence of additive noise occurs often in engineering. The parameters of 
interest are typically the tone's amplitude, phase, frequency and its time 
derivatives. Here our primary interest is in real-time estimation of dynamic 
doppler frequency, assuming constant amplitude. Estimation of the received 
phase is considered to be of secondary importance. 

Accurate frequency estimation and tracking is of fundamental concern in 
the design of Global Positioning System (GPS) receivers observing signals that 
exhibit high dynamics. GPS receivers process radio signals received from a 
constellation of 18 satellites in order to determine the receiver's position 
and velocity. Each signal is the sum of two L-band carriers, biphase 
modulated by pseudorandom codes and information streams. We assume that the 
modulation is known and has been removed from the signal. This assumption is 
valid when the ground station establishes a direct link with the GPS 
satellites, in parallel with the link to the transponder aboard the 
maneuvering target. Decoded data from the direct link can then be used to 
wipe data from the transponder link, effectively transforming these signals 
into unmodulated RF tones. The problem then reduces to estimating the 
parameters of two sinusoidal signals in additive noise. Here we concentrate 
on one of these tones, namely the LI carrier at a frequency of 1.575 GHz. 

In a previous paper [1], an estimator structure based on the maximum 
likelihood estimator of code delay and doppler frequency over a single symbol 
interval was analyzed. At signal-to-noise ratios above 30 dB-Hz, this 
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estimator performed well in the presence of considerable acceleration and 
jerk, even though no attempt was made to compensate for the effects of these 
trajectory components. The estimation error was close to the Cramer-Rao bound 
above 30 dB-Hz, below which rapid performance deterioration occurred, leading 
to loss-of-lock at about 28 dB-Hz. 

In the current application, the data have been removed from the carrier, 
hence we assume the estimator observes a sinusoid in the presence of additive 
noise. The frequency of the received sinusoid varies with time according to 
the trajectory of the vehicle being tracked. In principle, the estimator is 
free to observe the signal over several symbol intervals, thus decreasing both 
the estimation errors and the loss-of-lock threshold. However, rapid changes 
in the trajectory limit the maximum time interval during which the estimated 
parameters remain constant. Thus, a classical trade-off between noise 
suppression and doppler tracking results. 

Four different estimation techniques will be compared in the following 
sections. First we consider an extension of the "approximate maximum 
likelihood" technique described in [1]. This technique involves a significant 
extension of previous results, hence it will be developed in some detail. 
Second, an extended Kalman filter with parameters matched to the dynamics of 
the signal is described and analyzed, again in some detail due to the novelty 
of the high-dynamic application. The third technique employs a cross-product 
frequency control loop to obtain frequency estimates. Finally, the fourth 
approach examines a standard digital phase-locked loop. In all cases, 
estimator performance is compared on the basis of root mean squared estimation 
error and probability of loss-of-lock, while tracking a common simulated 
high-dynamic trajectory. 
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This article is subdivided into six sections. Section II describes the 


common trajectory used in the simulations. Section III is devoted to the 
development of the appropriate signal and noise models. The estimation 

algorithms are described and analyzed in Section IV, and the simulation 
results are compared in Section V. Summary and conclusions are presented in 
Section VI. Finally, the Cramer-Rao bound on the estimation error for the 
case of unknown acceleration is derived in the Appendix. 
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I I . COMMON TRAJECTORY 


The performance of the various estimators is evaluated on the basis of 
their ability to track a common trajectory. The trajectory chosen for 
simulation consists of positive- and negative-going jerk pulses of 0.5-sec 
duration and magnitude of 100 g/sec, separated by two seconds of constant 
acceleration, as shown in Fig. 1(a). The corresponding acceleration and 
velocity trajectories are shown in Figs. 1(b) and 1(c). The initial 
conditions for acceleration were chosen for symmetric 25-g excursions. The 
velocity trajectory is converted to an equivalent doppler frequency trajectory 
as 


f d ( t ) 



(Hz) 


( 1 ) 


which shows the one-to-one correspondence between the velocity v,(t) of the 
physical trajectory and the instantaneous frequency of the received signal. 
Here v^(t) is the doppler velocity, f^ denotes the carrier frequency, and 
c is the speed of light. Examination of Fig. 1(c) reveals that the doppler 
velocity, hence doppler frequency, of this trajectory can be well approximated 
by a piecewise linear model over time intervals on the order of 0.1 sec, even 
when a 100 g/sec jerk pulse is applied. This observation forms the basis for 
reducing the complexity of the simulated maximum likelihood estimator, as 
described in Section IV. 
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The common input to the simulation programs consists of samples of 
in-phase and quadrature sinusoids with phases proportional to the integral of 


the simulated frequency, plus independent additive noise samples, 
properties of these input samples are easily derived from the continuous 
of the received waveform. Therefore we begin by developing a suitable 
for the trajectory-modulated received signal. 


The 

model 

model 


6 



III. THE RECEIVED SIGNAL 


In the absence of noise, the received signal at the antenna may be 
represented as 

( _ _ j« c t) 

s(t)=Re|/2s(t)e j (2a) 

s(t) = Ae j0(t) (2b) 

where s(t) is the signal complex envelope, A is its amplitude, 

and where the phase process 0(t) is defined in terms of the doppler 

frequency process f d (t) as 


0(t) = 2ir f f d (0 dC 
—00 


(3) 


This trajectory-modulated signal is observed in the presence of a zero mean, 
stationary, narrow-band Gaussian noise process n(t), with representation 

I - ~ jw c t I 

n(t) = Re W2 n(t) e j (4) 

where n(t) is the complex envelope of the narrowband process. If the bandwidth 
of the noise process greatly exceeds that of the signal, then the covariance 
function of the complex envelope may be approximated as 


K~(t) A E[n(t + x) n* ( t ) ] ^ N n $(x) (5) 

where Nq is the two-sided spectral level of n(t). Henceforth, we assume 
that Nq is known. One can show that the two-sided spectral level of the 
corresponding real bandpass process is Nq/ 2. In the current application. 
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the carrier is first removed by mixing the received signal with sinusoids at 
the carrier frequency, yielding in-phase and quadrature components. 
Typically, both the in-phase and quadrature signals are first sampled at a 

high rate (mega-samples per second), and subsequently input to digital 

/ 

accumulators. The output of the accumulators, viewed as complex samples, can 
be modeled as 



(i+l)T fi 


L 


r(t ) dt = s . 

1 


+ n. 

i 


( 6 ) 


where r(t) = s(t) + n(t) is the complex envelope of the received signal. In 

the subsequent simulations, the effective integration time, T g , will be exactly 

2 msec in duration. The effects of spectral roll-off due to the effective 

integration will be ignored in this paper. If T is short compared to the 

s 

variation in s(t), then by the mean-value theorem there exists a t^ in the 

i-th T g interval such that = s(t^). For suitably small T g , we may regard 

t^ as the center of the i-th integration interval. 


Derivation of noise sample statistics is straightforward. Since n(t) is 

a complex zero mean Gaussian process, ru is a complex zero-mean random 
variable with covariance 


(i+1 )T S 


E<vt> - <V ' 2 


I -H/ 


(k+l)T fi 


N„ 


dt n(t ) n" ( t „ ) =i5. 
2 1 2 t l 


f~ ik 
x s 


(7) 
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is the 


Here the overbar denotes the expectation operator, and 
Kroenecker delta. Since all of the samples are pairwise uncorrelated, the 
joint density of any N distinct samples can be expressed as the product of the 
individual densities. This property is fundamental to much of our subsequent 
analysis. 
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IV. DESCRIPTION AND ANALYSIS 


In this section, the various frequency estimation algorithms are 
described and analyzed. Some of the algorithms also provide useful phase 
estimates; however, we shall not address the phase tracking problem 
explicitly. Analytical expressions for estimator performance are obtained 
where possible. In all cases, estimator performance is evaluated by means of 
numerical simulations. 


MAXIMUM LIKELIHOOD (ML) ESTIMATOR 


The maximum likelihood (ML) estimates of the signal parameters are those 
values that simultaneously maximize the conditional joint probability density 
of the observation vector, conditioned on the signal parameters. If the 
statistical distribution of the signal parameters within some uncertainty 
interval is not known, then maximum likelihood estimation yields the smallest 
estimation error variance. For each estimate, the observation vector consists 
of consecutive samples obtained over a time interval that is short compared to 
the characteristic timescale of the trajectory variations. We begin by 
developing a representation for the received phase process. 


Near a point tg in the domain of the trajectory function the frequency 
process can be represented in terms of a Taylor series expansion as 


f d ( t) 


-E 

0 k=0 


’ , (t - v* 

*k k! 


K 3t K 


t=t. 


(8a) 


(8b) 
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This representation is valid where all derivatives exist. For the trajectory 
function illustrated in Fig. 1, we can let f fe = 0 for k > 2, everywhere 
except at the points of discontinuity. As an example, near the point t^ = 0, 
representing the center of an observation interval. 


0(t) = 0 A + 2ir(f A t + 


(v 


JL 

2 




) 


(9) 


The coefficient f^ represents the doppler frequency in Hz, f^ its first 

derivative with units of Hz/sec, and f^ its second derivative with respect 
2 

to time (Hz/sec ). These coefficients correspond to the velocity, accelera- 
tion, and jerk of the physical trajectory. Henceforth we shall use the vector 
notation f = (f^, f^, f^) to denote these signal parameters. 


a. Maximum Likelihood Estimator Structure. 


The maximum likelihood estimator uses the conditional probability density 
of the observable vector, conditioned on the parameters of interest. Since 
the signal in Eq. (2b) is completely specified by the parameters A, 0 q, and f, 
it is equivalent to condition on the signal samples as 


p(fli) 



(N/2)-l 


E 




; N even 


( 10 ) 


The maximum likelihood estimates of the 
values that simultaneously maximize this 


signal parameters are those 
conditional joint density. 
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Equivalently, we can maximize any monotonically increasing functional of the 
conditional joint density, such as the logarithm. Taking the natural log of 
Eq. (10) and discarding terms which contain no information about the signal 
yields the "log-likelihood" function 



(ID 


Here r is an observation vector of length N, and r^, s^ are the samples 
defined in Eq. (6). With no loss in generality, the observation interval is 
assumed to be centered around the origin. Using Eq. (9) in Eq. (2b), we write 
the signal samples as 


s . 

i 



jgid) 


gi(f) = 2irf () t i 
h i (f 1 ,f 2 ) = 2ir 


+ h.(f l ,f 2 ) 



Substituting Eq. (12) into Eq. (11) yields 


(12a) 

(12b) 


(12c) 



(13) 


in terms of the modified samples 4 r^ exp(-jtu). This is the function to 
be maximized. Following the development in [1], we recall that for any 
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complex x, Re[exp(-j0g)x] is maximized with respect to 0g when 0g = arg(x), 
taking on the value |x|. Thus, 

e Q = arg 


z.e 

1 


-J 2 *Vi 


(14) 


Substituting into Eq. (13) yields 

- A 2 ^ (15) 

Since the amplitude is assumed to be unknown, its estimate A is found by 
differentiating Eq. (15) with respect to A, equating to zero and solving: 


max A(r) = ( rp I/ 2 A 

9 0 V 0 


-j2irf 0 ti 


z^e 


A 


1 

N 


E j 2irf 0 t i 

Z i e 


(16) 


Using this value in Eq. (15) yields 


max 

r o ,A 


*<■> ■ © E 


-j2irf 0 t. 


(17) 


Recalling that z. contains the desired parameters f^ and f 2 » it follows 

that the maximum likelihood estimates (fg, f^» ^ 2 ^ are those values that 
simultaneously maximize Eq. (17). Letting t^ = (i + (l/2))T g and dividing 
by the constant coefficient, we can equivalently maximize 


L(fg,f lt f 2 ) = 


(N/ 2 )— 1 


E »,■ 

i=-N/2 


_ -j 2irf o iT s 

e 


(18) 
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The function L(f Q , f-, f 2 ) is proportional to the squared magnitude of the 

Fourier transform of the modified sequence {5L}. This interpretation is 
useful for implementation, since efficient FFT algorithms can be employed to 
carry out the maximization. In general, estimator complexity depends on the 
number of significant terms in the Taylor series expansion of the received 
phase, which in turn depends on the trajectory. Thus, there is a direct link 
between the trajectory and the effective dimensions of the maximum likelihood 
estimator. 


b. Estimator Performance. 

In the current application, estimation of the received frequency is 

emphasized. The other parameters are regarded as "nuisance parameters" which 

may nevertheless have to be estimated to achieve best performance. The 

estimates are typically subject to threshold effects characteristic of 

nonlinear estimators, which means that below some critical carrier-to-noise 
2 

ratio (CNR A A /Nq) catastrophic deterioration in estimator performance 
begins to occur. Above threshold, estimator performance can be assessed by 
means of Cramer-Rao bounds. Below threshold, however, the breakdown mechanism 
must be modeled accurately to arrive at useful error estimates. We begin by 
developing some concepts that are necessary for understanding estimator 
performance. 

The maximum likelihood estimates of frequency and its time derivatives 
are those values f^, f^, and f^ that simultaneously maximize the real 
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function L(f). Let f'^, 2- = 0,1,2 denote the true parameter values, and let 
G(Af) denote the function 

(N/2)-l 

G(Af ) = N -1 ‘y ^ exp {-jg.(M)) (19) 

i=-N/2 

where Af^ = f^ - f£, and g^Af) is defined in Eq. (12b). We observe that 
2 

|G(Af)| is proportional to L(f) centered over f' when the variance of the 

2 

noise samples approaches zero. The projections of |G(Af)| onto cartesian 

coordinates defined by the signal parameters are shown in Figs. 2(a-c) where we 

let G(Afg) denote G(Afg,0,0), and so on. Note that the peak of the function 

2 

| G ( Af ^ ) | always occurs at the true parameter values, and that its shape 
depends only on the deviation from the true value: the coordinates of the peak 
track the temporal variations in the signal parameters, while the shape of the 
function remains unchanged. Therefore at high CNR, we expect the estimates to 
be near the peak along each coordinate, and the estimation errors to be 
related to the main lobe dimensions. 

2 

The effects of acceleration and jerk errors on | G(Af q ) | are 
shown in Figs. 2(d-f). Evidently, acceleration errors result in decreased 
peak amplitude (Fig. 2(d)), whereas jerk errors give rise to a shift in the 
peak as well (Fig. 2(e)). Significantly reduced and shifted peak values can 
result when both acceleration and jerk errors are present simultaneously, as 
shown in Fig. 2(f). 
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I G(A, 1>| 2 |g( A f Q ) | 
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Figure 2(c). [ G(Af 2 ) I 2 as a Function of Deviation from the 

True Derivative of Frequency Rate 







|G(Af 0 )| 2 |G< A f 0 i f 
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Figure 2(e). The Effect of Jerk Error on |G(Afg)|2 



Figure 2(f). The Effect of Both Acceleration and Jerk Error 
on |G(Af 0 )|2 
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The addition of noise impairs the estimator's ability to resolve the 
coordinates of the peak. The effects of additive noise can be quantified by 
means of the noise coherence function. For the case under consideration, we 
denote the effective three-dimensional transformed noise process by 
n^Cfg, f^, fg). Letting the signal component approach zero in Eq. (18), and 
using Eqs. (7) and (19), we obtain 


E{n 3 (f + Af)n*(f)} = £ £ «.< 

i & 


-j[g.(f+M)-g a (f)] 



G(Af) 


( 20 ) 


where 

) 

n 3 (f ) n. e 

i 

Thus, Eq. (20) exhibits the same dependence on Af as the function G(Af) 
defined in Eq. (19). 


An analytic expression for the projection of G(Af) along the frequency 
coordinate follows by letting t^= (i + (1/2)) T s in Eq. (19): 


G(Af Q ) = N 


sin (irAf-NT ) 
u s 

sin (irAf rt T ) 

U s 


( 21 ) 


This is recognized as the Fourier transform of an NT second pulse-train. 
However, for projections along the other coordinates, no closed form 
expressions could be found. 


The Cramer-Rao bound for estimating the frequency of a pure tone without 
phase information is well known. Using the results of [2], and converting to 


our notation, we have 
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var(4fo) 2 W(S) 


t\(n 2 - 1) 


; N > 2, known frequency rate 


( 22 ) 


This bound applies to unbiased estimators. The variance of estimation error 
is inversely proportional to CNR, and for large N, inversely proportional to 
the cube of the total observation time. If the frequency drifts with some 
acceleration, then Eq. (22) may still apply at high CNR, provided the 
acceleration does not change significantly from one observation interval to 
the next, and assuming that it has been estimated accurately in the past. If, 
however, the acceleration changes significantly over the timescale of an 
observation interval, so that it cannot be assumed known, then Eq. (22) does 
not apply. In that case, the Cramer-Rao bound for estimating frequency in the 
presence of an unknown frequency rate should be used. This bound, derived in 
the Appendix, has the form 


var 


Wfo> - 16 ( 2 ! 2 )(a°) tV ’ " 


> 2, unknown frequency rate 


(23) 


for large N. 


At low CNRs, estimator performance is critically linked to the total num- 
ber of "coherence cells" in the observation volume. Referring to Figs. 2(a-c), 
the "coherence-length" 6f^ along the ith coordinate is defined as the distance 
from the peak to the first minimum of the function |G(Af)| along that 
coordinate. We immediately conclude from Eq. (21) that Sf^ = l/T g N. The 
values of Sf^ and Sf^ can be computed accurately for any value of N. 
Thus we may think of associating a distinct random variable with each coherence 
cell, statistically independent of all others. 
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Although the statistical distribution of the signal parameters is not 
known, we assume that with high probability their true value is contained in 
some "uncertainty volume" whose center and dimensions can be specified. This 
uncertainty volume may be thought of as a region in with dimensions 
(Fq» F lf F2), such that with high probability 




for each d. We shall show that the estimator threshold tends to decrease with 
decreasing volume, provided the true parameters remain within this volume. 


The total number of independent random variables contained in an 
observation volume V = FgF]^ is 



( 24 ) 


Thus we may think of partitioning the observation volume into M disjoint 

coherence cells, one of which contains signal plus noise, while the remaining 

(M-l) contain only noise. The noise samples associated with different 

coherence cells are independent complex Gaussian random variables, each with 

2 

zero mean and variance o = N(N_/T ). The mean value due to the 

U s 

signal is AN. The probability density of the magnitude of signal plus noise 
is the Rician density function 



( 25 ) 
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while the noise sample magnitudes are distributed according to the Rayleigh 
density as 

, x lim , x 2x , 2. 2x /n ,. 

p 2' x ' = A ^, 0 P-j/x) = ~2 exp V-x /o ) (26) 

a 

In Eq. (25), Iq(») denotes the modified Bessel function of zero order. 

Since the maximum likelihood estimator selects the coordinates of the largest 
peak as its estimate, large errors can occur if the magnitude of a noise 
sample in one or more coherence cells exceeds the magnitude of signal plus 
noise. This event is called an "outlier." Outliers occur with probability 
q = 1 - p, where 

M-l 

dx (27) 

2 2 2 

Note that p and q depend only on the ratio (NA) /a = A T /N_ and M. The 

s U 

behavior of q as a function of CNR is shown in Fig. 3 for NT = 100 msec and 

s 

various M. It is apparent that for 100 msec observation times the probability 
of selecting the wrong observation volume is small for CNRs greater than 
23 dB-Hz, while at lower CNRs outliers are virtually certain. At any CNR, q 
is an increasing function of M, hence it is advantageous to keep M small in 
order to avoid frequent outliers, leading to large estimation errors. 

Given that an outlier occurred, the estimation error variance can be found 
by assuming that each point within the observation region has equal a priori 
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probability of being chosen, hence the frequency estimate may take on any 
value between -Fq/ 2 and Fq/ 2. Therefore, 


varCAfglcell error) = F” 1 



(Hz 2 ) 


(28) 


The high CNR result of Eq. (22) may be interpreted as the variance of the 
estimation error given that the correct coherence-cell was chosen, that is. 



Figure 3. Outlier Probability q as a Function of CNR 
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varCAfglno cell error). At any CNR, the variance of frequency estimation 
error is given by 

var(Afg) = q var(fg|cell error) + p varCf^Ino cell error) (29a) 


which can be bounded in terms of the underlying parameters as 


var (Af o ) > q 




(29b) 


The standard deviation of the estimation error as a function of CNR is 
shown in Fig. 4, for N = 50, T g = 2 msec and M = 160. For this example, the 
boundary between the high and low CNR regions is roughly 24 dB-Hz. This 
boundary need not correspond exactly to the loss-of-lock threshold, which we 
define as the CNR at which the estimator loses lock ten percent of the time. 



15 17 19 21 23 25 27 29 


CNR, dB-Hz 


Figure 4. RMS Frequency Error and its Components (ML) 




c. Simulation Results. 

Estimator performance was verified by means of numerical simulations. 
The received samples were obtained from the common trajectory by sampling at 
2-msec intervals, generating complex signal samples with the proper phase 
variation, and adding independent complex noise samples. In this 
implementation, N consecutive samples were processed by the estimator at one 
time, and the resulting parameter estimates referred to the center of the 
observation interval. These current estimates were used to predict the center 
of the observation volume for the following interval. 

A discrete approximation to the maximum likelihood estimator may be 
obtained by evaluating L(f) at uniformly spaced points along each coordinate 
of the uncertainty volume, with spacing fine enough to resolve the main lobe. 
Using this approach, efficient fast-Fourier transform (FFT) techniques may be 
employed to reduce the computational burden. In the current application, the 
simulated trajectory is dominated by 100-g/sec pulses of 0.5-sec duration. 
Since the sampling rate is fixed, the maximum change in acceleration increases 
with the total number of consecutive samples N. However, during constant 
acceleration, the rms estimation error decreases with increasing N. Thus we 
expect that for a given set of dynamic conditions an optimum value of N exists 
that provides high CNR without incurring excessive change in acceleration, 
which would require extending the acceleration uncertainty interval. 

The dimensions of the uncertainty volume and the number of discrete grid 

points along each coordinate are fundamental simulation parameters. Referring 

to Fig. 2(c) we note that for N = 50 and T =2 msec, the main lobe along the 

s 
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2 

Af^ coordinate extends to about +40000 Hz/sec . We also observe from Fig. 2(e) 

~ 2 

that an error in jerk of +100 g/sec = ±5000 Hz/sec causes only a slight 

2 

decrease in the peak of | G( Af Q ) | and a relatively small shift in frequency. 
Since we are primarily interested in estimator performance near threshold 

where frequency errors tend to be large, it seems reasonable to ignore the 
variation in jerk, and effectively assume it to be zero throughout the 

trajectory. Of course, this approximation impacts estimator performance at 
CNRs well above threshold, where the dominant source of frequency error now 
becomes the shift in the peak due to the non-zero jerk. 

Since we are projecting the estimates ahead in time to the center of the 

following observation interval, the maximum change in acceleration that can 

occur in NT = 100 msec due to ±100-g/sec jerk is +10 g = +500 Hz/sec. Thus 

s 

we let (F^/2) > 10 g, and partition F 1 into 2.5-g increments covering a 
range of +15 g, including a margin. The resulting 13-point grid samples 

G(Af^) at 125-Hz/sec intervals, easily resolving the main lobe. 

Finally, we recall that Fq = 1/T g = 500 Hz. Since the first zero 

along the frequency coordinate occurs at 10 Hz, the main lobe is well resolved 
by oversampling at roughly 2 Hz, implying the use of a 256-point FFT (the 
power of two nearest 250). Thus, we append 206 zeros to the 50 samples, and 
perform 256-point FFTs for each of 13 acceleration increments. From Eq. (24), 
using the value Sf^ = 700 Hz, the total number of independent random 

variables in this observation volume is M = [(500/10) + 1] [ (1500/700) + 1] = 
160, which is the value used in the computation of the outlier probabilities. 
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Representative sample sequences of frequency and frequency rate 
estimation errors are shown in Figs. 5(a) and 5(b), at a CNR of 23 dB-Hz. 
Note the highly quantized nature of the frequency rate estimation error, due 
to the 2.5-g quantization of the acceleration estimate. At such low 
carrier-to-noise ratios the estimator sometimes loses lock before completing 
the entire trajectory. 



TIME, sec 


Figure 5(a)* Instantaneous Frequency Error as a 
Function of Time (ML) 



TIME, sec 

Figure 5(b). Instantaneous Frequency Rate Error as a 
Function of Time (ML) 
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Figure 6(a). Estimated RMS Frequency Errors as 
Functions of CNR (ML) 



Figure 6(b). Estimated RMS Frequency Rate Errors as 
Functions of CNR (ML) 

Estimates of the rms frequency estimation errors are shown in Fig. 6(a) 
for N = 30, 50, and 80 using 2-g, 2.5-g, and 3~g acceleration quantization, 
respectively. The saturation in the estimation error at CNRs above 26 dB-Hz 
for N * 50 and 80 is attributed to the shift in the peak of the likelihood 
function introduced by the non-zero jerk pulses. These errors tend to 
increase with N because the magnitude of the shift is proportional to the 
observation interval. Simulation results for N = 50 were also obtained using 
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the true value of the jerk in order to bound the improvement in the frequency 
estimates that could be obtained if the jerk were estimated as well. The 
dotted curve in Fig. 6(a) presents the results. Apparently, no improvement 
is possible at or below threshold. Above threshold, however, frequency 
estimation errors could be reduced by the simultaneous estimation of jerk, 
acceleration, and velocity components. From Fig. 6(b) we note that the 
frequency rate estimates are not affected by mismodeling the jerk, at least in 
the range of the CNRs considered. 

If the predicted frequency for the upcoming observation volume is in 
error by more than +250 Hz, then the signal parameters most likely fall 
outside the succeeding observation volume, in which case loss-of-lock occurs. 
The loss-of-lock probability p is approximated by counting the number of times 
the estimator loses lock in 100 independent simulations. Estimates of the 
loss-of-lock probability are shown in Fig. 7 for N = 30, 50, and 80 samples 
per estimate. Note that loss-of-lock probabilities do not improve 
significantly if the true jerk is used, as shown by the dotted curve. 

Finally we observe that estimator performance could be improved somewhat 
by searching over a subset of the uncertainty volume. This approach is 
motivated by the observation that the maximum change in frequency is limited 
to +125 Hz, due to 25-g acceleration applied for 0.1 sec. Hence the total 
number of degrees of freedom becomes M * 26 x 3 = 78, which may result in a 
0.3-dB reduction in the threshold. However, such fine-tuning is possible only 
because the trajectory parameters are well known. If that were not the case, 
this minor improvement may well be cancelled by the increased design margin 
required to compensate for the greater initial uncertainties. 
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Figure 7. Estimated Loss-of-Lock Probabil- 
ities as Functions of CNR (ML) 


EXTENDED KALMAN FILTER (EKF ) 

Here we consider the application of an extended Kalman filter (EKF) to 
the problem of estimating the parameters of the received signal, driven by the 
common trajectory. A fundamental difference between this and the maximum 
likelihood approach is that the EKF provides instantaneous estimates after 
each new sample based on the latest sample value and previous estimates, 
whereas the maximum likelihood estimator provides estimates of average 
parameter values after processing a large number of samples. We begin by 
deriving the structure of the EKF estimator. 
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a. Extended Kalman Filter Estimator Structure. 


In this derivation we adopt vector notation, both for convenience and to 
conform to established practice. The received samples defined in Eq. (6) are 
now expressed as 


r(k) 


A sin [0(k ) ] 
A cos [0(k)] 


+ n(k) 


(30) 


where n (k) = [nj(k), n^(k)] is a zero-mean Gaussian vector, with the 
subscripts I and Q denoting in-phase and quadrature components, respectively. 
These correspond to the real and imaginary components of the complex noise 
sample defined in Eq. (6). From here on, the amplitude A is set to unity as 
it is assumed known and need not be estimated by the EKF. As before, we have 


E[n(k)] = 0 


(31a) 


E[n(k) n T a)] = 



k * ft 


(31b) 


2 

with I denoting the 2x2 identity matrix, and a = N n /2T . The 

n us 

extended Kalman filter operates on the phase 0(k) which is a sampled version 
of the integral of the frequency trajectory. The phase can be modeled as an 
n-th order polynomial whose derivatives constitute the components of the 
state-vector x(k) as follows: 


0(k) = !i T x(k) ; = [1, 0, 0] (32a) 

x(k + 1) = x(k) + v(k) (32b) 
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In the above, $ denotes the state transition matrix and v(k) is the disturbance 
vector. For a fourth-order EKF, x^(k) becomes 


x T (k) = [0(k), ugOO, “lOO* 


(33) 


where the various derivatives of 0(k) are denoted by w^Ck), to^(k), 
and {^(k) » respectively. Here each < 0 ^ is 2ir times as great as the 
corresponding f^ defined previously. The disturbance term models the 
random changes in the parameters due to dynamics. Since consecutive components 
of the state vector are integrals of the succeeding ones, repeated integration 
over the observation time yields 


« 2 ( k + 1) = + v 4^ k ^ 


(34a) 


»^(k + 1) = to 1 (k) + T g to 2 (k) + v 3 (k) 


(34b) 


a) Q (k + 1) = u) Q (k) + T g <o 1 (k) + co 2 (k) + v 2 (k) 


(34c) 


T 2 T 3 

0(k + 1) = 0(k) + Tw 0 (k) + -| o> 1 (k) + -| to 2 (k) + v^k) 


(34d) 


where 


kT 


'i‘ k) - / ¥(, > d¥ ! 


1 — 1, ..., 4 


(k-l)T 


(35) 
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and Y(t) denotes the fourth derivative of the continuous version of the 
phase. Assuming that Y(t) is a zero-mean white process with spectral level 
N y , we get 


E[v 2 (k)] = T = o 2 T 2 
4 2 s y s 


(36) 


2 

with denoting the variance of the sampled version of Y(t), and 

the level of the power spectral density of Y(t). The system matrix £ and 
the covariance matrix Q of the disturbance vector v(k) can be derived from 
Eqs. (34) and (35), as in [3]: 


£ = 


1 

0 

0 

0 


Q = N T 
3 y s 


T T 2 /2 T 3 /6 

s s s 

1 T T 2 /2 

s s 

0 1 T 

s 


0 0 

1 

J 


T 6 /252 

s 

T 5 /72 

s 

T 4 /30 

S 

T 3 / 24 
s 

T 5 /72 

T 4 /20 

3 

T J /8 

2 

T /6 

s 

S 

s 

s 

T A /30 

S 

T 3 /8 

s 

T 2 /3 

s 

t n 

s 

T 3 /24 

s 

T 2 /6 

s 

T /2 
s 

1 


(37a) 


(37b) 
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The following equations can then be derived for the estimation of the state 


vector x(k), as in [4]: 

x(k|k) = x(k | k - 1) + L(k)[r(k) - h(x(k|k - 1))] (38a) 

x(k + l|k) = ± x(k|k) (38b) 

Q(k) = H T (k) £(k|k - 1) H(k) + R (38c) 

L(k) = £(k | k - 1) H(k) Q _1 (k) (38d) 

E(k|k) = £(k|k - 1) - £(k|k - 1) H(k) 

X (H T (k) E(k |k - 1) H(k) + E)" 1 E T (k) £(k|k - 1) (38e) 

£(k + l|k) = a 2 £ £(k|k) + Q (38f ) 


where R = I, a is some weighting coefficient greater than one, 
T 

and H (k) is the linearization of the function h(x(k)), i,e., 


H T (k) = h(x) 


x = x(k | k-1 ) 


(39) 


where 


h(x | k) = 


sin(Q/^ x(k)) 
cos(fi;^ x(k)) 
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The constant a is used to adjust the weighting of past data in order to 
speed up filter response in transient situations. 

b. Estimator Performance* 

The performance of the fourth-order extended Kalman filter operating in 
dynamic environments and in the presence of additive noise can be bounded by 
means of the Cramer-Rao bounds. At high CNRs and mild dynamics the Cramer-Rao 
bound of Eq. (22) should apply, since this filter is an approximation to an 
optimum tracking algorithm. For high dynamics, where the variation in 
acceleration is significant, the bound for unknown acceleration, Eq. (23), 
should apply. Since this bound requires knowledge of the effective 
integration time, this quantity has to be determined. Response to dynamics, 
however, will only be assessed by means of the simulation results, except to 
note that dynamics tends to move the performance curves closer to the bound 
for unknown acceleration. 

c. Simulation Results. 

The performance of the EKF while tracking the common simulated trajectory 
has been evaluated. Both third- and fourth-order filters were considered. 
Since our primary interest is in evaluating tracking ability rather than 
acquisition, we assume that the initial conditions of the trajectory are 
precisely known. Thus the state components of the EKF are initially matched 
to the true trajectory. The exponential weighting coefficient a controls 
the memory of the filter. To reduce the effects of additive noise, a should 
be set to one in order to average over all previous data. In a dynamic 
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environment a is adjusted such that the present estimate is dominated by the 
past l/(a - 1) data points. For fixed sampling times, this establishes the 
equivalent integration time to be used in the Cramer-Rao bounds. An 
exponential weighting coefficient of 1.027 was used for the fourth-order EKF, 
as this value yielded the best performance. This corresponds to an effective 
memory of roughly 40 samples. 

A typical sample sequence of the instantaneous phase error is shown in 
Fig. 8(a) at a CNR of 26 dB-Hz, using a fourth-order EKF. It can be seen that 
the trajectory has little effect on the instantaneous phase error, while the 
instantaneous frequency and frequency rate errors tend to peak near the 
beginning and the end of each jerk pulse (Figs. 8 [b ] and [c]). This occurs 
because the peaks due to the trajectory are exaggerated by the presence of 
noise due to nonlinear effects within the estimator. Estimates of the 
loss-of-lock probability for phase and frequency estimates are shown in 
Figs. 9(a) and (b) for third- and fourth-order EKFs. Defining the threshold 
to be the CNR that yields 10% probability of loss-of-lock, it is clear that 
the fourth-order filter has a phase threshold of about 24.5 dB-Hz and a 
frequency threshold of roughly 24 dB-Hz. The rms phase, frequency, and 
frequency rate errors as functions of CNR are shown in Figs. 10(a), (b), and 
(c). At threshold, the minimum rms phase error is 0.4 rad, the minimum rms 
frequency error is about 3 Hz, while the minimum rms frequency rate error is 
roughly 75 Hz/sec. Note that the filter configuration with the lowest 
threshold does not achieve the lowest rms estimation errors. This is because 
the filter parameters were selected for minimizing loss-of-lock threshold 
rather than estimation error. 
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Figure 8(a). 


TIME, sec 

Instantaneous Phase Error as a Function of 
Time (EKF) 
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Figure 8(b). Instantaneous Frequency Error as a Function of 
Time (EKF ) 
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Figure 8(c). Instantaneous Acceleration Error as a Function 
of Time (EKF) 
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Figure 9(a). Estimated Phase Loss-of-Lock 

Probabilities as Functions of CNR (EKF) 


EKF PERFORMANCE 
O n = 4, a - 1.04 
□ n = 4, a = 1.03 
A n = 4, a = 1.027 
O n = 3, a = 1.06 
V n — 3, a = 1.05 

(n = ORDER) 
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Figure 9(b). Estimated Frequency Loss-of-Lock 

Probabilities as Functions of CNR (EKF) 
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CROSS-PRODUCT AUTOMATIC FREQUENCY CONTROL (CPAFC) LOOP 


In this approach a frequency discriminator is employed to track the 
received frequency. An automatic frequency control (AFC) configuration that 
employs a cross-product discriminator is also referred to as a 
"quadri-correlator" AFC loop. 


a. CPAFC Loop Estimator Structure. 


The discriminator output is represented in [5] by 


V k “ I k-1 Q k Q k-1 *k 


(40) 


where I and Q are the outputs of the in-phase and quadrature mixers. The 
performance of this loop operating in the presence of additive white Gaussian 
noise has been analyzed in [5] for a general loop filter. In our application 
the loop filter was chosen to be of the form 



(41) 


where z is the independent variable of the z-transform, = rd/T^, 

k n = rd^/T , d = 4B_ T/(r+l) f B T is the nominal bandwidth of the cross- 

product AFC loop, r is the damping parameter, and T the loop update time. 

s 

This loop is capable of tracking acceleration with zero steady-state frequency 
error. 
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b. 


Estimator Performance. 


In the absence of jerk and higher-order derivatives in the trajectory, 
the loop defined above tracks frequency variations with zero error in the 
absence of additive noise. However, jerk introduces a steady-state error that 
cannot be overcome with this loop. The addition of noise also degrades 
estimator performance, increasing the estimation error variance inversely with 
the carrier-to-noise ratio. Thus it makes sense to examine performance 
degradation due to either dynamics or additive noise separately, while the 
contribution of the other component is ignored. This allows us to assess loop 
performance when either component dominates. In general, both components 
contribute to the total estimation error. 


First we consider the steady-state frequency estimation error induced by 

a constant jerk of magnitude J. From [6], the steady-state frequency error 

3 

due to constant jerk of magnitude J (m/sec ) is 

(Hz) (42) 


s s 


= 5.25 - 
r 


r + 
4B 


which depends directly on the jerk, but inversely on the square of the nominal 
loop bandwidth. The variance of the frequency estimation error due to 
zero-mean additive Gaussian noise was found [6] to be 



(Hz 2 ) 


(43) 
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The error variance is seen to be inversely proportional to CNR and the square 
of the update time, but increases linearly with loop bandwidth. Comparing 
Eq. (42) with Eq. (43), we see that increasing the loop bandwidth yields 
improved steady-state performance at the cost of degraded noise performance. 
Thus, conflicting requirements must be balanced when both additive noise and 
severe trajectory variations are present. 

c. Simulation Results. 

The performance of the cross-product AFC loop was evaluated by means of 

numerical simulations. For the trajectory under consideration, it was found 

that a nominal loop bandwidth of = 8 Hz minimized the loss-of-lock 

threshold. At a carrier-to-noise ratio of 25 dB-Hz, with jerk of magnitude 

100 g/sec, this bandwidth yields a calculated steady-state frequency error 

of f = 22.6 Hz, and a noise-induced rms frequency error of roughly 
s s 

o £ = 18 Hz. The simulations confirmed these predictions. 

A typical frequency error trajectory is shown in Fig. 11 at a CNR of 
26 dB-Hz. The average frequency error appears to remain constant throughout 
the trajectory. Fig. 12(a) shows the estimated loss-of-lock performance of 
the loop for various bandwidths and damping ratios, while Fig. 12(b) displays 
the corresponding rms frequency errors as a function of CNR. The lowest 
loss-of-lock threshold appears to be 24.5 dB-Hz, achieved by a 7-Hz loop with 
a damping ratio of 2. The corresponding rms frequency error is roughly 40 Hz 
for this case. We observe that small variations in bandwidth do not affect 
the loss-of-lock performance of the AFC loop significantly, but seem to have a 
more serious effect on the rms estimation error. 
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TIME, sec 


Figure 11. Instantaneous Frequency Error as a Function of 
Time (CPAFC) 
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PHASE LOCKED LOOP (PLL) FREQUENCY ESTIMATOR 


Frequency estimates can also be obtained from the traditional digital 
phase locked loop. The block diagram of this loop is basically the same as 
that of the AFC loop, with the discriminator and loop filter functions 
redefined. Because of the dynamics inherent in the common trajectory, a 
type-III loop was simulated since it can track linear frequency variations 
with zero steady-state error. 

a. Phase Locked Loop Estimator Structure. 

The loop filter was chosen to be an "impulse invariant transformat ion" 
filter (IIT) due to its simplicity and desirable performance characteristics. 
The filter transfer function [6] is 

G 2 G 3 

FU) = ♦ rj ♦ (44) 

1 - z (1 - z ) 

where G, = rd/T , G„ = rd^/T , G_ = rkd^/T , and d = 4B t T (r - k)/r(r - k + 1). 
1 s 2 s o s Ls 

Here is again the nominal loop bandwidth, used to define the parameter d. 
Note that a second-order loop filter and the third-order numerically controlled 
oscillator used in the simulation [6] yield a fifth-order digital loop. 
However, this loop behaves like a type-III analog loop. For the simulations, 
a sinusoidal phase detector characteristic was assumed. Frequency estimates 
were obtained from the phase samples via the difference equation 
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(Hz) 


(45) 


f = 


0-0 - 
n n-1 


n 2irT 


where 0^ denotes the phase estimate sample at time nT g < 


b. Estimator Performance. 


For the parameter values defined above, we have from [6] that the 

3 

steady-state phase error due to a constant jerk of magnitude J (m/sec ) 


is 



co 

c 

c 


J n T 3 
0 s 

rk 


r(r - k + 1) 

4B T (r - k) 
L s 


(rad) 


(46) 


while the phase error variance due to noise is 

B l (rad 2 ) (47) 

The loop bandwidth should be chosen to minimize the sum of the estimation 
errors due to both trajectory and additive noise. 



c. Simulation Results. 


It was found by simulation that a nominal loop bandwidth of 43 Hz (r = 3, 

k ss 0.5) minimized the probability of losing lock for our trajectory. For 

this bandwidth with J = 100 g/sec, Eqs. (51) and (52) yield 0 gs = 0.46 rad and 

a = 0.095 rad at a CNR of 26 dB-Hz. The simulations compare favorably 
0 

with theory for two extreme cases, namely the case without thermal noise and 
the case without jerk. Estimates of the rms phase estimation error and the 
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frequency estimation error are shown in Figs, 13(a) and (b), respectively. 


Both estimation errors increase with decreasing CNR without much evidence of 


thresholding. However, we note from Fig. 13(c) that the probability of 


loss-of-lock exhibits a well-defined threshold around 26 dB-Hz. 
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Figure 13(a). Estimated RMS Phase Errors as 
Functions of CNR (PLL) 





PROBABILITY OF LOSS-OF-LOCK 



Figure 13(b). Estimated RMS Frequency Errors as 
Functions of CNR (PLL) 
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Figure 13(c). Estimated Loss-of-Lock Probabilities 
as Functions of CNR (PLL) 
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COMPARISON OF RESULTS 


As mentioned in the introduction, our primary interest is to estimate the 
frequency of the incoming signal. Other signal parameters, such as phase or 
frequency rate are considered to be of secondary importance here. Some of the 
algorithms naturally provide estimates of secondary parameters in the process 
of estimating the primary ones. For example, the maximum likelihood algorithm 
estimates the rate-of-change of frequency, the phase locked loop derives 
frequency estimates from phase estimates, while the extended Kalman filter 
estimates all three parameters simultaneously. In the following paragraphs we 
compare the various algorithms on the basis of their ability to track a common 
frequency trajectory, while commenting on their ability to estimate the other 
parameters. 

The loss-of-lock probability for each algorithm as a function of CNR is 
shown in Fig. 14(a). The parameter set that achieved the lowest threshold was 
selected in each case. Threshold is defined as the CNR at which the 
loss-of-lock probability for frequency estimation is 0.1. The approximate 
maximum likelihood estimator achieved the lowest threshold among the 
algorithms tested, namely 23 dB-Hz. The threshold for the extended Kalman 
filter was 24 dB-Hz, whereas the cross-product AFC loop and the phase-locked 
loop attained thresholds of 24.7 dB-Hz and 25.7 dB-Hz, respectively. The 
difference in thresholds is due to unequal sensitivity to severe dynamics 
among the various implementations. 

The root-mean-squared (rms) frequency estimation errors are displayed in 
Fig. 14(b), again as functions of CNR. For each algorithm, the parameter set 
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RMS FREQUENCY ERROR, Hz PROBABILITY OF LOSS-OF-LOCK 






corresponding to minimum threshold was used. Note the sharp increase in rms 
estimation error for the ML algorithm near the loss-of-lock threshold. The rms 
estimation errors of the algorithms differ at threshold, which is typically 
the lowest CNR of interest. The EKF achieves roughly 3-Hz rms error at 
24 dB-Hz, the PLL reaches 14-Hz rms at 25.7 dB-Hz, the ML yields 7 -Hz rms at 
23 dB-Hz , while the AFC algorithm operates with 40-Hz rms error at 
24.7 dB-Hz. These values represent the smallest estimation errors each 
algorithm can achieve at threshold, below which loss-of-lock will frequently 
occur. 

These algorithms should also be compared on the basis of operational 
comlexity, which we define as the average number of computations required per 
sample. The ML algorithm is the most complex, requiring the rapid computation 
of thirteen 256-point complex FFTs. The EKF requires approximately one order 
of magnitude fewer computations than the ML algorithm. The AFC and PLL 
implementations are the least complex, with the PLL algorithm requiring the 
fewest computations . 

The relevant estimator attributes are summarized in Table I. The last 
column indicates increasing relative complexity, ranked from one to four. The 
other columns display threshold, rms frequency estimation error at 26 dB-Hz 
(well above threshold for most estimators), and the availability of phase and 
frequency rate estimates. For many applications, the choice of a suitable 
estimator algorithm can be based on these attributes. 
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Table I. Overall Comparison of the Four Frequency Estimators 


Algorithm 

Threshold 

(dB-Hz) 

RMS Error (Hz) 

23 dB-Hz 26 dB-Hz 

Phase 

Estimate 

Frequency 

Rate 

Estimate 

Complexity 

ML 

23.0 

7 

1 

No 

Yes 

4 

EKF 

23.9 

3.5 

2.2 

Yes 

Yes 

3 

CPAFC 

24.7 

60 

31 

No 

Yes 

2 

PLL 

25.7 

20 

12 

Yes 

No 

1 
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VI . SUMMARY AND CONCLUSIONS 


In this paper, four different frequency estimation algorithms were 
presented and discussed. An approximate maximum likelihood algorithm, an 
extended Kalman filter, a cross-product AFC loop, and a phase locked loop were 
compared on the basis of tracking a common frequency trajectory that exhibits 
severe dynamics. The results are applicable to high-dynamic trajectories in 
general. The maximum likelihood approach was found to attain the lowest 
loss-of-lock threshold (23 dB-Hz), and also the lowest rms estimation errors 
above threshold. Although the performance of the extended Kalman filter was 
somewhat worse in both respects, it was able to operate with lower frequency 
estimation errors near threshold. The digital phase-locked loop performed 
well above threshold, but could not maintain lock reliably below about 
26 dB-Hz. The threshold for the cross-product AFC loop was somewhat lower 
than for the phase-locked loop, but its estimation errors above threshold were 
the greatest of the four algorithms tested. Based on these simulations, we 
conclude that frequency estimates approaching the theoretical bound for known 
accelerations can be achieved in dynamic environments comparable to our 
simulated trajectory by the approximate ML and the EKF algorithms. 
Performance in terms of loss-of-lock probability and estimation error will 
generally depend on the severity of the dynamics encountered, and on the 
extent to which the system parameters can be matched to the temporal 
characteristics of the dynamic trajectory. 
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APPENDIX 


DERIVATION OF CRAMER-RAO BOUNDS 


The Cramer-Rao (CR) bounds on estimation error variance are derived for the 
case of a sinusoid of unknown frequency, frequency rate, and phase. We assume 
that the following in-phase and quadrature measurements are available: 

r x (k) = A sin(e Q + + B Q t£) + n^k) 4 s^k) + n^k) 

r 2 (k) = A cos(9 0 + 03 Q t k + (3^) + n 2 (k) 4 s 2 (k) + n 2 (k) (Al) 

t^ = t^ + kT = (n^ + k) T ; k = 0, 1, ... N - 1 

where n^(k), n 2 (k) are components of the noise vector satisfying Eq. (31b). 

We wish to obtain CR bounds for the unbiased estimation of the unknown 

but constant parameters 9^, co^, and (3^. Denoting the parameter vector 
T 

(|3q» Qq) by a, and the conditional probability density function of the 

observation (conditioned on a) by f(R;a), R(k) = (r(k), 0 k <_ N - 1}, r(k) = 
T 

[r^(k), r 2 (k)] , and following the development in [2], the ij-th element of the 
Fisher information matrix J becomes 


J. . 

iJ 


= E 


!aS7 f ^> 



3s^(k) 3s^(k) 

3a . 3a . 

i J 


3s 2 (k) 3s 2 (k)| 

+ 3a . 3a . ( 

i J ) 


where a^ = |3 q, a 2 = Wq, a^ = 9 q. After some straightforward computation, we 
obtain 
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N-l 


n k=0 


n 

L 3 

n 

t 2 

n 


n 

L 2 

n 

t 

n 


n 


n 


(A2) 


Letting 


J »ii J - 

A 


b 

c 

d 


(A3) 


the elements of the matrix may be expressed in terms of N and n^ as 


a = T 4 {n 4 N + 4n 3 P + 6n 2 Q + 4n Q R + S} 


b = T 3 {n 3 N + 3n 2 P + 3n Q Q + R} 


c = T 2 {n 2 N + 2n Q P + Q} 


(A4 ) 


d = 


T {n Q N + P} 


e = P 


where 
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N-l 

P = k = N(N - 1 ) / 2 
k=0 


N-l 


Q = ^ k 2 = N(N - 1)(2N - 1 ) / 6 
k=0 


R 


=X > 3 


= N 2 (N - l) 2 /4 


(A5 ) 


N-l 

k=0 


k = N(N - 1 ) (2N - 1)(3N - 3N - l)/30 


The bound on the error variance of the estimate of a. 

i 

2 2 

(a /A ) multiplied by the i-th diagonal element of 
can be expressed as 


is then given by 
J„ . This matrix 


'N 




(ce - d ) 


(ae - c ) 


(ac - b ) 


with 

|J N I = a(ce - d 2 ) - b(be - cd) + c(bd - c 2 ) 


(A6 ) 
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Substituting Eqs. (A4) and (A5) into Eq. (A6) yields 


ce - d = CjT 


ae - c 2 = ( z m 2 C 1 + 4n Q C 2 + C 3 > T 4 


(A7 ) 


ac - b - ( n gC^ + ^ n g^2 + n Q^4 + 2n 0^5 + ^6^ ^ 


The explicit expression for |j^| and various £s can be obtained in 
terms of N as follows: 


| jr I = T 6 N 3 (N - l) 2 (N 4 + 2N 3 - 3N 2 - 8N - 4)/2160 (A8) 

C x = N 2 (N 2 - 1)/12 

C 2 = N 2 (N - l) 2 (N + 1)/12 

C 3 = N 2 (N - 1) ( 2N - 1) (8N 2 - 3N - 11)/180 

(A9 ) 

C 4 = N 2 (N - 1) (7N 3 - 23N 2 - 18N - 8)/60 

C 5 = N 2 (N - l) 2 (2N 3 - 3N 2 - 3N + 2)/120 

C = N 2 (N - l) 2 (8 IN 4 - 162N 3 + 105N 2 - 56N - 8)/240 

o 


The variance of the frequency estimate f^ = coq/2it, the parameter of greatest 
interest, then has the following bound 


var(f Q ) > (v 2 n 0 + + v Q ) 


A 2 T 2 A(2ir) 2 


(A10) 
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where 


*2 = 4C 1’ V 1 = 4C 2’ V 0 = C 3’ A = I J nI /T 


For the special case = 0, 


var(f-) > 


v o a 


0 - A 2 T 2 A(2*) 2 


(All) 


If N is so great that the highest-degree terms in N dominate in the 
expressions for and A, then 


, / N c 

var(f-) > -r 


96 


0 "\a 2 / T 3 N 3 (2ir) 2 


(A12) 


Again for n^ = 0, the error bound for [3 q is 


var(|3_) > 


N 2 (N 2 - 


0' - 12A 


n (£) 


(A13) 


which may be approximated for large N as 


\ ~ 180 a 
var(|3_ ) > 


0 2 4 5 

ATT 


(A14) 


2 

Noting that 3^ = 2w^, with co^ the frequency derivative in rad/sec , and also 
2 2 

that (a / A T ) = (N n /2P T ) where (P /N_) is the received signal power-to-noise 
s U c s c u 

spectral density ratio, we obtain 

var(i i> >^(r) <A15) 


62 



2 

in units of (Hz/sec) . Comparison with the bound in Eq. (23) 
performance bound for estimating f^ is a factor of 16 
presence of an unknown frequency rate, than for the known (or 
rate) case. 


shows that the 
greater in the 
zero frequency 
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